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Abstract 

The hydrostatic primitive equations (HPE) form the basis of most numerical 
weather, climate and global ocean circulation models. Analytical (not statistical) 
methods are used to find a scaling proportional to (Nu Ra Re) 1 ^ for the range of 
horizontal spatial sizes in HPE solutions, which is much broader than currently 
achievable computationally. The range of scales for the HPE is determined from 
an analytical bound on the time-averaged enstrophy of the horizontal circulation. 
This bound allows the formation of very small spatial scales, whose existence would 
excite unphysically large linear oscillation frequencies and gravity wave speeds. 

The hydrostatic primitive equations (HPE) have been the foundation of most numerical 
weather, climate and global ocean circulation calculations for many decades [H [2j [3J, IU 
In practice, modern computational power can handle integrations of these on global 
horizontal grids ranging in size between 15km and 60km, which correspond respectively 
to one-eighth degree and one-half degree in latitude and longitude at the equator. This 
limitation raises the long-standing question, "Can numerical simulations at these grid 
sizes adequately predict climate and other natural phenomena that occur on the much 
wider range of scales observed in Nature?" See Figure [TJ 




Fi gure 1: A NASA image [5] illustrates the large range of fluid scales that exist in atmospheric circu- 
lation. The oceanic range of scales is similar, but is not so easily observed. 



J D Gibbon and D D Holm (13th/ 09/ 2011) 



Enstrophy Scale Range of HPE 



2 



Important as it may be, this long-standing question is not addressed here. Rather, 
two questions are addressed associated with the HPE model itself. Namely, "What range 
of scales is available for solutions of the HPE?" and "What scaling law governs the size 
of horizontal HPE excitations in terms of the system parameters?" These dimensionless 
parameters are Nu, Ra and Re, associated with the names of Nusselt, Rayleigh and 
Reynolds, respectively. 

The HPE differ from the three-dimensional Navier-Stokes equations in incorporating 
rotation, stratification, and imposing vertical hydrostatic balance. The latter is often 
regarded as the most accurate of the various assumptions used in large-scale computations 
of the climate, weather and ocean circulation. The hydrostatic assumption determines 
the pressure from the weight of the water above a given point, independently of its state 
of motion. This changes the nature of the dynamics, because the vertical velocity is 
determined from incompressibility, rather than from its own evolution equation. 

Unlike the Navier-Stokes equations, solutions of the HPE have been proved to be 
regular by Cao and Titi [7]. Moreover, the HPE have also been shown to possess a 
global attractor [Sj. Although its solutions are regular, the HPE system may potentially 
possess a vast range of sizes of excitations [9j. While Kolmogorov introduced the Re 3 ^ 
scaling law for the range of spatial sizes of excitations in incompressible fluid flows by 
using statistical methods |10| . the present note will use analytical methods to show that 
a scaling law exists, proportional to (Nu Ra Re) 1 ^, for the range of horizontal spatial 
sizes in solutions of the HPE, with similar boundary conditions to those of Cao and Titi 
[7J. This result demonstrates that HPE excitations are possible at scales that are many 
orders of magnitude smaller than are possible in present numerical resolutions. 

A dimensionless version of the HPE may be expressed in terms of two sets of velocity 
vectors involving the horizontal velocities u, v and the vertical velocity w [11] 

V(x,y,z,t) = (u, v, ew), v = (u, v, 0). (1) 

Under the constraint of incompressibility, div V = 0, these satisfy 

e (d t + V ■ V) u - v = eRe^Au - d x P, (2) 

e (d t + V ■ V) v + u = eRe^Av - d y P. (3) 

Here e is the Rossby number, Re = UqL/v is the Reynolds number and P the pressure. 

As mentioned earlier, HPE has no evolution equation for the vertical velocity com- 
ponent w. Instead, this variable is determined (diagnosed) from the incompressibility 
condition, divV = 0. The z-derivative of the pressure field P and the dimensionless 
temperature enter through the equation for hydrostatic balance 

a 6 + d z P = . (4) 

The coefficient ao = (ea~ 1 a~ 2 )RaRe~ 2 arises from non-dimensionalization of the equa- 
tions. Here o = p/k is the Prandtl number (the ratio of viscosity v and thermal diffu- 
sivity k), Ra is the Rayleigh number, defined by Ra = gaToH 3 ^/^)' 1 , g is acceleration 
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of gravity, a is volumetric expansion coefficient, To is a typical temperature difference 
and a a = H/L is the aspect ratio of the cylindrical domain. When ([2]), ^ and Q are 
combined, an evolution equation for the hydrostatic velocity field v = {u, v, 0) results as 

e (d t + V ■ V) v + k x v + a kB = eRe~ 1 A v - VP, (5) 

which is taken in tandem with the incompressibility condition divV = 0. The dimen- 
sionless temperature G (the source of buoyancy) evolves according to 

(d t + V ■V)e = (aRe)- 1 AQ + q, (6) 

in which nondimensional q specifies heat sources, or sinks. The domain Q is taken to be 
a cylinder of radius L and height H. The vertical velocity and vertical flux of horizontal 
momentum both vanish on its flat upper and lower cylinder surfaces (z = 0, H). That is, 
w = and u z = v z = on the boundary. The variables are all taken to be periodic on 
the sides of the cylinder. 

Linearizing the HPE in ^ and ([6]), and their non-hydrostatic equivalent (which has 
the dynamics of w restored), leads to well-known dispersion relations [13], which are 
illustrated in Figure |2j The essence of these dispersion curves is that without the frequency 
cut-off enforced by the buoyancy terms in the non-hydrostatic equations, the HPE admit 
unphysically high gravity wave frequencies at small scales. Moreover, these HPE gravity 
waves propagate at a fixed phase speed in the limit of small scales, while in reality gravity 
waves at these scales cease to propagate at all. 



1 An estimate for the resolution length 

Taking the inner product of the divergence-free velocity V with the motion equation ^ 
gives an equation for the rate of change of the kinetic energy of horizontal motion 

~ / \v\ 2 dV= [ (Re- l V ■ Av - a wQ) dV (7) 
2 dt Jq Jq 

in which dV is the volume element and surface terms integrate to zero under the present 
boundary conditions. For the Navier-Stokes equations it is normal practice to use the 
energy dissipation rate v (f n |<-o>| 2 <iV) based on the full vorticity uj = curlV to define 
a length scale called the Kolmogorov length |12j . The quantity J n \u}\ 2 dV is called the 
enstrophy and the angle brackets ( • ) denote the time average 

1 ' T 



(•)= lim - / (-)dt. (8) 







However, it is more appropriate in the hydrostatic approximation to use three-dimensional 
C = curl v and base a horizontal length scale on ( J n \ £ \ 2 dV) , since the vertical velocity w is 
diagnosed from the horizontal velocity dynamics. To determine this horizontal length scale 
from the evolution of the horizontal kinetic energy in ([7]), let us examine the Laplacian 
term 

r V ■ AvdV = - I u> ■ CdV (9) 
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Horizontal wavenumber kH 



Fi gure 2: This comparison of linear mode dispersion relations for the hydrostatic primitive equations 
(solid curves) with those of the exact nonhydrostatic equations (dashed curves) for oceanic conditions 
shows that the primitive equations admit very high fluctuation frequencies, especially at high horizontal 
wave numbers. In contrast, the dispersion relation for the nonhydrostatic equations limits properly to the 
buoyancy frequency, regardless of how high the horizontal wave number becomes. Oceanic parameters 
are taken as c s = 1500ms -1 , H = 10 3 5 m, iVo = 0.01s -1 where the appropriate normalizing length scale 
H is the mean ocean depth. The multiple curves correspond to different choices of vertical wave number, 
mH £ [0, 1, 3, 9, 27], increasing from the left. The value m = is the barotropic mode and the others are 
baroclinic. 



where the surface terms again vanish for our choice of boundary conditions. Note that £ 
is fully three dimensional, but its horizontal components vanish at the top and bottom of 
the cylinder. Two more integrations by parts give 

f u>-CdV= [ (|C| 2 + (div^) 2 ) dV> [ \C\ 2 dV (10) 
Jn Jn Jn 

and thus Q may equivalently be re-written as 

/ \v\ 2 dV <-Re- l f (CI 2 dV - a I wB dV . (11) 
Jn Jn Jn 



Id 
2dt 



Upon defining the vertical Nusselt number Nu as 

Nu:=-( [ wBdV), (12) 



n 



the time average of (11) may be written as 



(13) 



since the horizontal kinetic energy term vanishes in the limit as T — > oo. This bound on 
the time-averaged enstrophy of the horizontal circulation (J* |£| 2 <iV) yields a horizontal 



J D Gibbon and D D Holm (13th/ 09/ 2011) 



Enstrophy Scale Range of HPE 



5 



resolution length scale which emerges upon switching back into dimensional variables. 
Let C,dim be the dimensional version of C, ; that is, Cdim = L^UqC, for a typical horizontal 
velocity scale U . Then a resolution scale X Tes may be defined using the same approach 
as that used to find an analytical estimate of the inverse Kolmogorov scale for the Navier- 
Stokes equations. 



r4 X -4 
u /\ res 



-2 r-3 



-2 



\C\ 2 alV 



Re 1 



\C\ 2 dv 



(14) 



Thus, the main result obtained from ( 13 ) and ( 14 ) is an estimate for the range of horizontal 
scales, defined by the ratio LX^q S , as 



£A~e S < {ea' l a~ 2 NuRaRe) 



1/4 



(15) 



This bound incorporates all physical processes in their nondimensional forms. Estimated 
from the time-averaged enstrophy of the horizontal circulation, the ratio LAfes °f the 
domain size to the resolution scale provides an upper bound for the range of horizontal (not 
vertical) length scales. The hydrostatic approximation holds regardless of the magnitude 
of this ratio. 



Conclusion 



It is now time to put some numbers into the estimate in (15). For example, in regional 
flows in the ocean of depth H pa 10 5 km, aspect ratio a a = 10~ 2 , Prandtl number o pa 10 
and Rossby number e = 10~ 2 , one has £cr _1 a~ 2 ~ 10 1 . Thus, the range of scales (15) in 
this case may be written as 

L ^res ^ ( 10 NuRa Re) l/A . (16) 

The Rayleigh, Prandtl and Nusselt numbers usually appear in Rayleigh-Benard convec- 
tion in which Nu is observed to scale with Ra such that Nu ~ Ra 13 with variations around 
/3 = 1/3 : see [H] for a discussion of the state of the art for heat transfer and large scale 
dynamics in turbulent Rayleigh-Benard convection. However, the hydrostatic approxima- 
tion excludes deep convective processes, in which case Nu ~ 1 [15] . The Rayleigh-Benard 
/3-scaling for Nu would apply only at small vertical turbulence scales where the hydro- 
static approximation would be invalid. An important issue in oceanic simulations is to 
differentiate between mass flux and heat flux. Numerical simulations of ocean circulation 
must typically be corrected to prevent over-estimating the heat flux [16]. The need for 
this correction is another indication that the Nusselt number tends to be small in oceanic 
flows. 
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The sizes of Ra and Re for typical flows in the ocean are very large, when based 
on regional domain size and molecular values of viscosity and diffusivity of heat. For 
example, with H pa 10 3 ' 5 m and Nu ~ 1 

Ra = gaT^H^^vK)^ 1 

^ 10 1 10- 4 10°10 10 - B (10 6 10 7 ) « 10 20 - 5 (17) 

and Re = U H/{ua a ) pa 10 _1 10 3 - 5 (10 6 10 2 ) « 10 10 - 5 . According to these estimates, 
RaRe~ 2 = 0(1), and the coefficient a = {ea~ l a~ 2 )RaRe^ 2 ~ 10 1 ; so the range of 
scales is bounded by about eight orders of magnitude. That is, in this case, LA^es ~ 10 8 - 
This means that for a domain size of 400km at a depth of about 4km, the horizontal 



excitation scales could be as small as a few millimeters. In particular, the estimate (15) 
with Nu pa 1 and Ra ~ Re 2 yields 

L\^ s <(e*- l a- 2 ) l/4 ReV\ (18) 

which is close to the Kolmogorov range of scales in 3D. The very high linear wave fre- 
quencies associated with such small horizontal scales would preclude both the physical 
relevance and the computability of the HPE. The conclusion is that improving the reso- 
lution of HPE numerical solutions may tend to make their results less accurate and much 
more expensive to perform, because the nonlinear tendency toward much smaller spatial 
scales produces wave excitations of rapidly increasing linear frequency (as in Fig. 2) that 
would require reducing the time-step beyond the present limits of computability. Appar- 
ently, this fact is already recognized in practice, since the HPE are generally applied to 
climate simulations, but not to regional simulations. What this paper shows and empha- 
sizes is that unphysically small spatial scales can potentially be generated in HPE with 
molecular values for transport coefficients. In fact, modulo appropriate adaptations, the 
same range of scales would be found to hold for the nonhydrostatic equations, although 
we do not discuss it here because no proof of existence is available for them. 

Of course, numerical simulations of large-scale circulations in the ocean and atmo- 
sphere do not use the molecular values of viscosity and diffusivity. Instead, they introduce 
effective values for these quantities due to unresolved scales, associated with turbulent 'ed- 
dies'. These effective values are chosen essentially to make the Reynolds number at the 
horizontal grid scale Re(Ax) equal to unity. If the scaling Ra ~ Re 2 persists for these 
simulations and the Nusselt number at the grid scale is of order unity, then the numerical 
procedure of setting Re(Ax) = 1 might tend to properly resolve the hydrostatic excita- 
tions of the HPE. However, it may also be good practice in numerical simulations using 
the HPE to evaluate the dimensionless numbers at the vertical grid scale Nu(Az) and 
Ra(Az) corresponding to the other physical aspects of the HPE. Further study of the 
scaling law Ra ~ Re 2 for various regimes of ocean and atmosphere circulation might also 
be fruitful in determining local values of the ranges of scales. 
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